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Covariance Dynamics and Entanglement in Translation Invariant 
Linear Quantum Stochastic Networks 
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Abstract 

This paper is concerned with a translation invariant network of identical quantum stochastic systems subjected to external 
quantum noise. Each node of the network is directly coupled to a finite number of its neighbours. This network is modelled 
as an open quantum harmonic oscillator and is governed by a set of linear quantum stochastic differential equations (QSDEs). 

The dynamic variables of the network satisfy the canonical commutation relations (CCRs). Similar large-scale networks can 
be found, for example, in quantum metamaterials and optical lattices. Using spatial Fourier transform techniques, we obtain a 
sufficient condition for stability of the network in the case of finite interaction range, and consider a mean square performance 
index for the stable network in the thermodynamic limit. The Peres-Horodecki-Simon separability criterion is employed in 
order to obtain sufficient and necessary conditions for quantum entanglement of bipartite systems of nodes of the network in 
the Gaussian invariant state. The results on stability and entanglement are extended to the infinite chain of the linear quantum 
systems by letting the number of nodes go to infinity. A numerical example is provided to illustrate the results. 

I. INTRODUCTION 

Recently, there has been an increased interest in the study of dynamics of large-scale quantum networks. In particular, one 
can refer to the studies on artificial optical media known as quantum metamaterials; see, for example, [20], [14], [32], [19], 
[32], [30], These networks are composed of complex unit cells and are effectively homogeneous on the scale of relevant 
wavelengths, for example, in the microwave range. In comparison with natural solids, whose optical behaviour is specified 
by the quantum energy level configurations of the constituent atoms or molecules, the electromagnetic response of quantum 
metamaterials stems from the controllable resonant characteristics of their building elements such as the Josephson devices 
or optical cavities [20], [14], [19]. Quantum metamaterials are considered to be a promising approach to the implementation 
of quantum computer elements, such as quantum bits (qubits) and qubit registers, which can maintain quantum coherence 
over many cycles of their internal evolution [19], [29]. 

The machinery of linear quantum stochastic differential equations (QSDEs) provides a framework for modelling and 
analysis of a wide range of open quantum systems [9], [15], including those which arise in quantum optics. In the QSDEs, 
the quantum noise from the surroundings is modelled as a heat bath of external fields acting on a boson Fock space [15]. The 
class of linear QSDEs represents the Heisenberg evolution of pairs of conjugate operators in a multi-mode open quantum 
harmonic oscillator which is coupled to external bosonic fields. In this framework, the modelling of translation invariant 
networks with field coupling between the subsystems and stability analysis for certain classes of nonlinear open quantum 
systems have been addressed, for example, in [27], [18], [23]. 

For modelling purposes, the analysis of homogeneous large-scale networks can be simplified by using periodic boundary 
conditions (PBCs) [19], [27]. The PBCs rely on negligibility of boundary effects in a fragment of a translation invariant 
network consisting of a sufficiently large number of subsystems. This approximation technique is used for lattice models of 
interacting particle systems in statistical physics, including, for example, the Ising model of ferromagnetism [12]. 

Quantum decoherence, as the loss of quantum nature of the system state towards classical states, is described by using 
entanglement measures and correlations between the nodes of bipartite systems [22]. Among purely quantum mechanical 
phenomena, entanglement (or inseparability) plays a crucial role in the fields of quantum information and quantum 
computation [13], [8], Various measures have been presented for entanglement; see, for example, [8]. A criterion for checking 
whether a given system state is separable was proposed in [17]. This condition was shown to be necessary and sufficient 
for separability of Gaussian states in two-mode quantum harmonic oscillators [24], 

In this paper, we consider large fragments of a translation invariant network of identical linear quantum stochastic systems 
endowed with the PBCs. The nodes of the network are directly coupled to each other (through bilateral energy transfer) within 
a finite interaction range. This network is modelled by linear QSDEs based on Hamiltonian and coupling parametrization of 
the corresponding multi-mode open quantum harmonic oscillator. The dynamic variables of the network satisfy the canonical 
commutation relations (CCRs). Using spatial Fourier transforms [27] and linear matrix inequalities (LMIs), we present 
sufficient conditions for stability of the network, thus extending the previously considered nearest neighbour interaction 
setting to the case when each node of the network is coupled to a finite number of its neighbours. We also consider a 
mean square performance index for the stable network in the thermodynamic limit. In the case when the network consists of 
one-mode oscillators, we obtain sufficient and necessary conditions for quantum entanglement between nodes of the network 
in the Gaussian invariant state, based on the Peres-Horodecki-Simon separability criteria. By letting the number of nodes go 
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to infinity, the results on stability and entanglement are extended to the infinite network of linear quantum stochastic systems. 
For the infinite network in the Gaussian invariant state, it is shown that the entanglement disappears for nodes which are 
sufficiently far away from each other. However, we provide a numerical example of a large sample of long fragments of a 
stable network with randomly generated parameters, where the bipartite entanglement between nodes is present within the 
interaction range. 

The rest of the paper is organized as follows. Section |H| provides principal notation. Section [HI] specifies the class of open 
quantum systems being considered. Section |IV] describes a model for a fragment of a translation invariant network of such 
systems with PBCs. Section [V] employs spatial Fourier transforms in order to obtain a more tractable description for the 
covariance dynamics of the network. Section VI discusses a stability condition and a mean square performance functional. 
Section [VH| p resents conditions for entanglement between nodes of the translation invariant network of one-mode oscillators. 
Section [VHI| extends the results on stability and entanglement to the infinite chain of such systems. Section |Ix| presents a 
numerical example to illustrate the results. Section |X| makes concluding remarks. Appendices [A][B] and [C] provide subsidiary 
material on a closed-form solution of the Lyapunov equation and entanglement criterion of bipartite systems, and some facts 
about matrices, respectively. 


II. NOTATION 

Unless specified otherwise, vectors are organized as columns, and the transpose (-) T acts on matrices with operator-valued 
entries as if the latter were scalars. For a vector X of operators X\ ,..., X r and a vector Y of operators Y \,..., Y s , the commutator 
matrix is defined as an (r x s)-matrix [X,T t ] :=XY t — (TA t ) t whose (j,k)th entry is the commutator [X n Y^] := XjY^ — Y^Xj 
of the operators Xj and Y^. Also, (■) '' := ((•) # ) T denotes the transpose of the entry-wise operator adjoint (■)*. In application 
to complex matrices, (•)' reduces to the complex conjugate transpose (•)* := ((-)) T . Furthermore, § r , A r and H,-: = 8 r + iA r 
denote the subspaces of real symmetric, real antisymmetric and complex Hermitian matrices of order r, respectively, with 
i := \f—\ the imaginary unit. Also, I, denotes the identity matrix of order r, positive (semi-) definiteness of matrices is 
denoted by Op) A, and (g) is the tensor product of spaces or operators (for example, the Kronecker product of matrices). 
The generalized (Moore-Penrose) inverse of a matrix M is denoted by M + , and r (M) denotes the spectral radius of M. The 
Kronecker delta is denoted by <5^, and D := {z £ C : |z| = 1} is the unit circle in the complex plain. The floor function is 
denoted by [-J. Also, Eq := Tr(p|) denotes the quantum expectation of a quantum variable q (or a matrix of such variables) 
over a density operator p which specifies the underlying quantum state. For matrices of quantum variables, the expectation 
is evaluated entry-wise. 


III. LINEAR QUANTUM STOCHASTIC SYSTEMS 

We consider a quantum stochastic system interacting with external boson fields. The system has N subsystems with 
associated ^-dimensional vectors Xq, ... ,X^~i of dynamic variables which satisfy the CCRs 


[Y,X t ] = 2/0, 


0 := In ® ©, 



( 1 ) 


Here, 0 is a 
0=4® L°! 


block diagonal joint CCR matrix, where 0 £ A 2 „ is a nonsingular matrix which is usually of the form 
q . The system variables evolve in time according to the QSDE 


dX=(i[H,X]-^BJB T &- l X^dt+BdW, 



Wv-iJ 


( 2 ) 


Here, Wo,...,Wv-i are 2;;;-dimensional vectors of quantum Wiener processes with a positive semi-definite ltd matrix £2 £ 

IHL m- 

dWjdW? = S jk Odt, -Q := hm + U, (3) 


. The matrix B £ R 2nNx2mN i n (2 


is 


where J := /,„ £ 

the system-field coupling operators on the system variables by 


related to a matrix M £ M 2mNx2nN 0 f linear dependence of 


B = 2 0M t . (4) 

The term \BJB J 0 l X in the drift of the QSDE (2ji is the Gorini-Kossakowski-Sudarshan-Lindblad decoherence 
superoperator [5], [10] which acts on the system variables and is associated with the system-field interaction. Also, H 
is the Hamiltonian which describes the self-energy of the system and is usually represented as a function of the system 
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variables. In the case of an open quantum harmonic oscillator [3], [4], the Hamiltonian H is a quadratic function of the 
system variables 

1 1 N ~ x 

H := -X t RX = £ xjR jk X k , ( 5 ) 

1 1 j,k =0 

where the matrix R := {R jk)o^j,k<N G §2 „n is formed from (2 n x 2/z)-blocks Rjg. By substituting (|5j» into (|2j) and using the 
CCRs 0. it follows that the QSDE takes the form of a linear QSDE 

dX = AXdt + BdW, (6) 

where the system matrix A £ m 2nN><2 nN j s gj ven by 

A:=20R-^BJB T 0~ l . ( 7 ) 

IV. LINEAR QUANTUM STOCHASTIC NETWORK WITH PERIODIC BOUNDARY CONDITIONS 

Suppose the open quantum harmonic oscillator of the previous section is a fragment of a translation invariant network 
which is organised as a one-dimensional chain of identical linear quantum stochastic systems numbered by k = 0,... ,N— 1. 
Each node of the network interacts with the corresponding external boson field and hence, the joint network-field coupling 
matrix M in 0 is block diagonal: 

M = I n ®M, (8) 

where M £ R mx ". Furthermore, the nodes are directly coupled to each other within a finite interaction range d. The fragment 
of the chain is assumed to be large enough in the sense that N > 2d, and is endowed with the PBCs, thus having a ring 
topology. A particular case of nearest neighbour interaction is depicted in Fig. |T| In the case of an arbitrary interaction range 



Fig. 1. A finite fragment of a translation invariant network of open quantum systems with direct coupling between the nearest neighbours. Also shown 
are the quantum noises which drive the nodes of the network. 


d^l, the Hamiltonian H in (Jsj) is completely specified by matrices Rg = R 


= R J P £ 


r>2nx2n 


, with £ = 0, ± 1,..., ±d, as 


TV-1 


H :=^ £ (xj £ ). 

z ]=0 v i=-d ' 


(9) 


Here, j — £ is computed modulo N in accordance with the PBCs, and hence, the corresponding matrix R is block circulant. 
The matrix Ro £ § 2 ,, specifies the free Hamiltonian for each node, while R± s describe the energy coupling between the nodes 
which are at a distance s = 1 ,d from each other (with the more distant nodes in the network not being directly coupled). 
In view of the block diagonal structure of the matrices 0 and M in 0 and 0, it follows from 0,0 and 0, that the 
QSDE 0 is representable as a set of coupled QSDEs for the dynamic variables of the nodes of the network: 


d Xj = £ A e X mod(j _ w dt + BdWj, j = ().... ,N — 1. 

e=-d 

Here, similarly to R, the network dynamics matrix A in 0 also has a block circulant structure, with 

_ J 2&R 0 -\BJB t ®~ 1 if £ = 0 
^ \ 2 ®R e if e ± 0 ’ 

B :=2©M t . 


( 10 ) 


( 11 ) 

( 12 ) 
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V. SPATIAL FOURIER TRANSFORMS AND COVARIANCE DYNAMICS OF THE NETWORK 
The dynamics ( [T()| of the nodes of the network can be studied in the spatial frequency domain. Similarly to [27], we will 
use the discrete Fourier transforms (DFTs) of the quantum processes and 144 over the spatial subscript k = 0. ...,1V — 1: 


Jfc=0 

N-l 


m) ■■= E z~ k w k ( t ), 


k =0 


where z £ Uv and 


f zniK 

Un := {e N : k = 0,...,N—l} 

is the set of /Vth roots of unity. The quantum process 3£ z (t) is evolved in time by the QSDE 


(13) 

(14) 

(15) 


N -1 

d = E z ~ JdX J 

y'=o 

JV-l / d 

= E Z M E 'd^^mod(y-t ; A')dt + BdWy 

7 =o v=-<* 

= t <^Cdr + /id'#[ ; 


where the matrix A z is defined in terms of the matrices GD 


as 


:= E *> 


(16) 


(17) 


k=—d 


and B is given by (121. Here, use is made of (13 i and ( | 1 4[ > together with the PBCs and the well-known properties of DFTs 
due to the assumption that z £ U#. Furthermore, it follows from the CCRs G that 

St'v) = E z~ ivk l x j,Xk] = 210 E (v/z)‘ = 2 i8 zv N® 

j,k =0 k=0 

for all z,v £ Ujv. By a similar reasoning, (|3]» implies that 

dWAWj =N5 zv £h\t. (18) 

We will need the following matrix of second-order cross-moments of the quantum processes SE Z and 3L V given by S 

y z ,v(t) :=E(^T-(t)^;(0 t ). (19) 

Here and in what follows, the quantum expectation E( ) is over the tensor product p := CJ 0 U of the initial state GT of the 
network and the vacuum state v of the external fields. 

Lemma 1: For any Nth roots of unity z,v£Uaj from (15», the matrix ,5P z . v (t) in (19i satisfies the differential Lyapunov 
equation 

= .<zi z ,9’ zy + + N8 ZV BL1B J (20) 

for any t ^ 0, where the matrix sf is given by Furthermore, suppose sz/ z is Hurwitz for all z, £ U,v • Then, for any 
Z,v £ Uyv, there exists a limit 

J^vH:= ; hm J?z,v(t) (21) 

which is a unique solution of the algebraic Lyapunov equation 

+NS zv BQB t = 0. (22) 

Proof: A combination of the QSDE ( | 1 6j i with the quantum Ito formula leads to 

d^, v = dE(^ v t ) 

= E(d 3f z ^ + + dJQd.^) 

= + ^z,v< +N5- r BDB T )dr, 

which implies the ODE (20 1 . Here, use is made of ( p~8] > and the fact that the forward increments of the quantum Wiener 
process in the vacuum state are uncorrelated with the adapted processes, and hence, K(,JJ~ z c\"//f,') = 0. The solution of (201 
is described by 

(23) 


y zv (t) =e , ^ v (0)e , <+AS n , f 

Jo 


e T - <;/ -/LQ/i T e T < dr 
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for any t ^ 0. Therefore, under the assumption that sf z is Hurwitz for all z £ U,y, the solution (23 i has the following limit 


r+°° 


y,Joo)=N8- r I e r ^BQ.B J e Ts ^dr 
Jo 

which is a unique steady-state solution of ( p0| reducing to ( |22} . 

Under the assumptions of Lemma [T| the solution of ( [22} is representable as 

'¥’zA°°)=Nd zv y z , 

for all z,v £ U,v. where y z = yf 0 is a unique solution of the algebraic Lyapunov equation 

£/,y,+y z y* + bub' = o. 


(24) 


(25) 


The function U 9 z >->• S z is the spatial spectral density [27] which encodes the covariance structure of the dynamic variables 
of the network in the invariant Gaussian quantum state (in the limit of infinite time and infinite network size). Also note 
that the cross-moments of the dynamic variables at the /th and kth nodes of the network can be recovered from the matrices 
© by applying the inverse DFT to ([13}: 


E(xj(t)x k (t) r ) = T E £ z j v~ k y z (t)y v {ty 


z,veUjv 


4 E 


(26) 


Z,VGUa( 


for all 0 ^ j,k < N. In fact, the right-hand side of (26i is the two-dimensional inverse DFT of the matrices ^,1/v 

VI. STABILITY AND MEAN SQUARE PERFORMANCE OF THE NETWORK 

Recall that the network matrix A in 0. which corresponds to the set of QSDEs ( [TO}, has a block circulant structure 
specified by the matrices HD- Therefore, A is Hurwitz if and only if the matrix sf z in (| 1 7} is Hurwitz for all z £ ILJ.y. The 
latter property holds for any fragment length N ^ 1 if and only if 


maxrle^) < 1 

zeU 


(27) 


(since s4 z is a continuous function of z over the unit circle, and the set U,v=i Ujy is dense in U). In what follows, the 
translation invariant network, described in Section IV is called stable if it satisfies (27 i. The following lemma provides a 
sufficient condition for stability of the network under consideration. 

Lemma 2: Suppose there exist positive definite matrices y, Q £ § 2 ,, satisfying the LMI 


A 0 y+yA^+Q y a 

& -hn 0 

A T 0 — jjhnd 


=4 0 , 


where A £ . 


2n/ And j s an aU xili ar y matrix defined in terms of (11 1 by 

A := [A- d ... A_i Ai ... A d \. 


Then the quantum network being considered is stable in the sense of ( [27} . 
Proof: In view of the Schur complement lemma [31] and ([29}, the LMI 


implies that 


A 0 y+yAl+Q+y 2 + 2 dY,A j A T j 40 . 

yo 


(28) 


(29) 


(30) 


For any z £ U, the vector A := [z d ,...,z,z 1 ,...,z d ] £ satisfies |A| 2 ^ 2d, and hence, AA* =<! 2dhd- Therefore, 

application of the completing-the-square technique leads to 


£ A jZ ~ j y +^ £ A] z j 42dY,AjA T j+y 2 , 

yo yo yo 


(31) 


where use is made of the relations Y*yoAjZ J = A(A®/ 2 h) and A(A(& h,f){A(A® h n ))* = A((AA*) (gi/ 2 „)A T 2dAA T . A 
combination of (30 1 and ( |3 1 } implies that 

A 0 y+yAg + £ Aj Z -jy+y £ A] z j + Q =* o, 


yo 


yo 


5 










which is equivalent to 


^/ z y+y^+Q40. 


(32) 


Since 39 and Q in (28 i were assumed to be positive definite, (32 1 implies that is Hurwitz for all z £ U, which proves 

■ 

Note that the condition for Aq to be Hurwitz is a prerequisite for feasibility of the LMI ( |28j ). Now, assuming that the 
stability condition ( [27] ) is satisfied, we will consider the following mean square performance measure [27] for the network 
of size N at time t ^ 0: 


At-1 


^v(0:=- E E Xftfoj-M). 


(33) 


j,k =0 


Here, o k is a given 


2 n x 2n 


valued sequence which satisfies = oj for all integers k and specifies a real symmetric block 


Toeplitz weighting matrix (<Jj~k)os^j,k<N- The block Toeplitz structure of the weighting matrix in (33 i corresponds to the 
translation invariance of the quantum network being considered. A matrix-valued map U 9 z H E z = £*, defined by 


E ; := £ z ~ k a k , 

k=-o o 


(34) 


describes the spectral density of the weighting sequence. In order to ensure the absolute convergence of the series in (34 1 , it 
is assumed that ||o*|| < +°°, which also makes E- a continuous function of z. The fulfillment of the condition E z 7= 0 

for all z £ U is necessary and sufficient for (Oj -k)o<j.k<N 7 0 to hold for all /V 7 1; see, for example [ 6 ], In this case, the 
sum on the right-hand side of ( |33] > is a positive semi-definite operator, and hence, <4%- 7 0 . 

Lemma 3: Suppose the stability condition ( [27] is satisfied. Then for any given length N of the network fragment, the cost 
functional in (331 has the following infinite time horizon limit 


1 


<&(“):= lim = E Tr^z)^). 


zGU^v 


Here, 


S*>=- I (-£) 

k=\-N 


z k Ok, 


(35) 


(36) 


and the matrix -9’ z is the unique solution of the algebraic Lyupunov equation (25 i. 

Proof: The relation ( |26| allows the expectation on the right-hand side of ( [33] to be represented in the spatial frequency 
domain as 


E(Xj(t) T a^ k X k (t)) = Tr(cjJ_,E(X ; (0^(0 T )) 


72 E zly 


N 2 


(37) 


z,vgUa) 


where use is also made of the property 07 1 = a k of the weighting sequence. Hence, in view of (21 1 , (24 1 and (36 1 , the 
infinite time horizon limit of Sfa takes the form 


4H = z3 E E z>v 


A 3 


j,k =0 
N -1 


jy 2 E E ^ kr ^ r (°k-j'9 z ) 

j,k =0 zeUiv 

1 


= - E Tr(E N (z)y z ), 


zZ U,v 


which establishes 

VII. BIPARTITE ENTANGLEMENT CRITERION FOR THE NETWORK 
For any fixed but otherwise arbitrary nodes j ^ k of the translation invariant network under consideration, we denote by 


t (f) — m 
X k (t) 


(38) 


the augmented vector of dynamic variable of the corresponding bipartite system. In a stable network driven by external 
fields in the vacuum state, the mean values of the dynamic variables asymptotically vanish and will be considered to be zero 
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in what follows. Since the joint unitary evolution of the network and its environment preserves the CCRs (jT|, the quantum 
covariance matrix of the vector T in ([38]) takes the form 


S(t) :=E(Y(f)Y(f) T ) = R eS{t)+i 


0 

0 


0 

0 


Suit) 

Suit) 


Sn{t) 
$22 (0 


>= 0, 


(39) 


where the off-diagonal blocks are real matrices satisfying S], = So i . Note that in view of (261, the matrix S is representable 
in terms of the matrices ([T9| whose evolution is described by Lemma |T| 


S(t) 


1 

Jf2 


E 

z,vGUw 


(z/v)J 

z k v~j 


z j v~ k ' 

{z/v) k 


®y zy {t). 


(40) 


Gaussian quantum states [16] are completely specified by the covariance matrix of the quantum variables and their mean 
vector, with the latter not affecting whether the state is separable or entangled. Therefore, separability or entanglement of 
such states can be studied in terms of the covariance matrices [28]. Furthermore, the Gaussian nature of system states is 
inherited by subsystems and is preserved by linear QSDEs driven by vacuum boson fields (provided the initial state is 
Gaussian). For what follows, we assume that the network consists of one-mode oscillators (that is, n = 1) and is initialized 
at a Gaussian state. Then by applying the Peres-Horodecki-Simon criterion [24], it follows that the above described bipartite 
system (formed from the /th and kth nodes of the network) is in a separable Gaussian state if and only if 


A(f) :=ReS(f) + i 


0 

0 


0 

-0 


—S [t ) + i 


’0 0 


\Sn(t) 

S 12 W 

0 -20 


bi 2 (t) T 

Soif) 


0, 


(41) 


where S(f) is the quantum covariance (4 x 4)-matrix in ( |39| . An equivalent alternative form of the separability criterion is 
provided by Lemma[4]of Appendix | b| which reduces to checking the sign of detA(t) instead of verifying the (4x4)-matrix 
inequality in ( [41) . More precisely, according to Lemma [4j the bipartite system is in an entangled Gaussian state if and only 
if the quantum covariance (2 x 2)-matrices .S) ] and S 22 for the constituent nodes j and k are positive definite and 


det(5 2 2-5l 2 5 n 1 5 1 2)<0. 


This is closely related to the following identity (see part (iii) of Appendix |C| for details): 

det A(f) = detSn det(S 22 - 3 , } 2 5' I H |.S) 2 ) 

= detS 22 det(Sn — 5 i 2 5 22 + 5’) 2 )- 


Theorem 1: Suppose the translation invariant network consists of N one-mode open quantum harmonic oscillators and is 
stable in the sense of ( [27) . Then the invariant Gaussian state of the bipartite system at the nodes j f^k is separable if and 
only if the infinite time horizon limit of the matrix A (?) in (41 1 satisfies 


A(°°) := lim A(f) = V 


y z z>- k y z 






)= 0 , 


(42) 


where .56. is the unique solution of the algebraic Lyapunov equation (25 1 . 

Proof: In view of Lemma [T] the fulfillment of the stability condition ( |27) implies (24 1 whose combination with (40 1 
leads to 


5(°°) := lim Sit) = — V 


. 7 - z j ~ k y z 


The latter allows the separability criterion (41 1 to be represented for the corresponding limit matrix A(°°) in the form d42). 


A closed-form solution of the algebraic Lyapunov equation ( |25) for the one-mode case (which is considered in Theorem[Tj) 
is presented in Appendix [A] 
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Fig. 2. An infinite chain of directly coupled linear open quantum systems with nearest neighbour interaction. Also shown are the quantum noises (numbered 
from left to right). 


VIII. INFINITE CHAIN OF LINEAR QUANTUM STOCHASTIC SYSTEMS 


We will now consider the infinite one-dimensional chain of directly coupled linear quantum stochastic systems whose 
finite fragments were introduced in Section IV A particular case of nearest neighbour interaction is depicted in Fig. [2] 
The results of Sections VI and |VII| which employ the representation of the translation invariant network in the spatial 
frequency domain (see Section [V[), can be extended to the infinite network case by letting the fragment length N go to 
infinity. The mean square cost functional S’n, discussed in Section VI leads to a performance index for the infinite network. 
More precisely, similarly to [27, Theorem 3], under the stability condition (27 1 , the infinite time horizon value <£/v(°°) in 
(|35|> has the following limit as the fragment length N tends to infinity: 


lim 

TV—>+°o 


M°°) = ^=r.<f Tr(X^) 

2m J u 


dz 

z 


(43) 


Here, E, is the spectral density ( [34] ) of the absolutely summable weighting sequence. The resulting cost functional in ( |43j ) 
corresponds to the thermodynamic limit in equilibrium statistical mechanics [21], The results of Theorem |T| on separability 
of states for bipartite subsystems (in the network of one-mode oscillators) can be extended to the infinite chain as follows. 

Theorem 2: Suppose the infinite chain of one-mode open quantum harmonic oscillators satisfies the stability condition 
(|27[. Then the invariant Gaussian state of the bipartite system at the nodes j ^ k is separable if and only if 




z j - k . r A 

Jze u 

_z k - j y z 

.9z 


A ‘ 2ni , 

where SP Z is the unique solution of the algebraic Lyapunov equation 
Proof: Under the stability condition \21\, the solution -9' z of 


dz 

— f 0, 

z 


(44) 


is analytic (and hence, continuous) with respect to 


Z over a neighbourhood of the unit circle U. Hence, the matrix Aat(°°) : = A(°°) in (42 1 , which is the Riemann sum of an 
appropriate integral, converges to the matrix A in (44 1 as N —» +°°. ■ 

Under the assumptions of Theorem [ 2 J the matrix A in (44 1 is representable as 


A = 


So 


(45) 


S U 

in terms of the Fourier coefficients of the spatial spectral density Z/' 7 \ 

Sf := i z r ~ 1 <9 z dz. 

2ni J ze u 

The latter describe the appropriate cross-covariance matrices for the invariant Gaussian state of the network and, in view of 
the Riemann-Lebesgue lemma, satisfy lim/ >00 .SV = 0. Therefore, if So A 0 (and hence. So >~ 0, see Appendix [c]), then (45) 
implies that A ^ 0 for all j — k large enough. This means separability (that is, the absence of entanglement) for all nodes 
j and k which are sufficiently distant from each other. 


IX. ILLUSTRATIVE EXAMPLE 


We have generated a sample of 100 sets of random coefficients of the QSDE ( fT0[ > for a translation invariant network of 
N = 400 one-mode (n = 1) oscillators subject to the stability condition of Lemma IT] The interaction range in the network 
was d = 8, and the dimension of the driving quantum noises was 2m = 2. For each of these samples, detA(°°) from (42 1 


and the log-negativity index [26] were calculated for the bipartite systems at nodes j of the network. The sample means 
and ranges of values of these entanglement measures versus a:= j — k are depicted in Fig. [3] The mean values show that 
it is quite probable for the nodes of the network to be entangled within the interaction range (that is, for |a| ^ 8 in this 
example). According to Fig. [3] the entanglement vanishes for sufficiently distant nodes (with «| >8 in this case). 
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Fig. 3. Mean values (solid lines) and ranges of values (colour filled domains) for detA(°°) and the log-negativity versus a := j — k for nodes j and k 
in 100 randomly generated stable translation invariant networks of 400 one-mode oscillators. Negative values of detA(°°) and the corresponding positive 
values of the log-negativity indicate the quantum entanglement between the nodes. 


X. CONCLUSION 

We have presented the modelling of translation invariant networks of directly coupled linear quantum stochastic systems 
with finite-range interaction. Using the spatial Fourier transforms and LMIs, we have provided conditions of stability for 
the network, and considered a mean square performance functional with block Toeplitz weights for the stable network. For 
networks of one-mode oscillators in a Gaussian state, we represented the Peres-Horodecki-Simon entanglement criterion 
(for bipartite systems at arbitrary nodes) in the spatial frequency domain. We have discussed an extension of the results on 
stability and entanglement to infinite chains of linear quantum systems by letting the network size go to infinity. We have 
also provided a numerical example which revealed the presence of entanglement between the nodes of the network within 
the interaction range. 
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Appendix 


A. Solution of the Algebraic Lyapunov Equation 
The solution SP Z of the algebraic Lyapunov equation 
of matrices [11] as 

vec(oS^) = —rvec (BQ.B r ), 

where the dependence of T on z is omitted for brevity, and 
sf z . Since srf z = \/ z for any z £ U in view of (17 1 , it follows from (A1 1 that 5A Z is a rational function of z which is analytic 

in a neighbourhood of the unit circle U, provided the stability condition (27 1 is satisfied. Indeed, the latter guarantees that 
the matrix srf z (£>£/ z is also Hurwitz (and hence, nonsingular) for any z € U. Similarly to the approach in [1], the matrix F 
can be computed as 


can be computed by using the column-wise vectorization vec(-) 

r — K®^.)- 1 , (Ai) 

> := h„ <8> sX v + stf 7 ®h.n is the Kronecker sum of stf z and 


F=(I 2n ®D x )(hn®£?? n 


2n—2 


+ £ 2 <S>< 2 "" 3 - ... - (—1) 2 "£2»-i ®/2»), 

(A2) 

D \=srf z n - c\srf z n ~ x + c 2 srf z n ~ 2 - ... + c 2n I 2n , 

(A3) 

El :=&/ z +C\l2m 

(A4) 

Ek .=srf z Ek- 1 + Ck^2n i k = 2,..., 2 n 1 

(A5) 


where c; are the coefficients of the characteristic polynomial of that is 


det(A/ 2 „ — s/f) = X Zn +c\X~ n 1 +... + C 2 „. 

In the case when the nodes of the translation invariant network are one-mode quantum harmonic oscillators (that is, n = 1), 
the representation ( |A2| >-( |A5j > reduces to 

r = (I 2 ®D~ l )(I 2 <gi srf z — E\ ®/ 2 ), 

= / 2 (8)(D^V.)-£ 1 (g)(D- 1 ), (A6) 


D = srff + srffFrsrf z + h det g/ z , 
E\= srf z — LTr.eA. 


(A7) 


Furthermore, in this case, the Caley-Hamilton theorem leads to .<//} = .fi/,Tr.&/, — l 2 det///-. and hence, the matrix D in (A7 1 
takes the form 

D = Tr(^- + &/f)£/ z + (det^- — det^/ z )I 2 . 

In order to obtain a closed-form expression for ,jA, the matrix I) 1 can be computed as 


D 


1 

det I) 


(I 2 TiD-D). 


Here, 


det/5 = Tr(^A + ■sXfj dzt(X z I 2 — &/ z ) 

= Tr(j4 + £/*) 2 (A z -Tr(^)A z + det(^)), 
= 4(Rev z ) 2 (X z — v z X z + u z ), 


where v z := Tr.fA, u z := dct.xz), and X- := !™ z . Therefore, 


D = 


Re t ; :((vz ~ 21m u z )I 2 - sX-f) 
(Imw-) 2 — 2v z Rev z lmu z + 2(Rev z ) 2 u z 


This implies that 


y z = (^/ z bfib t - bfib t eJ ), 


where use is made of ( |A1[ ), ( |A6| > and matricization techniques. 

B. An Equivalent Bipartite Entanglement Criterion for Two-mode Quantum Harmonic Oscillators in Gaussian States 
Lemma 4: Suppose a two-mode quantum harmonic oscillator with the CCR matrix I 2 0 © is in a Gaussian state with the 
quantum covariance (4 x 4)-matrix S in |39] l. Then this state is entangled if and only if the matrix A in (41 1 satisfies 

det A < 0. (Bl) 

Proof: Positive semi-definiteness of a Hermitian matrix is equivalent to nonnegativeness of all its principal minors [2, 
p. 421], In view of (39i and ( |4T| >, all the principal minors of A up to order 3 inherit nonnegativeness from those of the 
quantum covariance matrix S 0 (see part (iii) of Appendix [C] for details). Therefore, violation of the property A ® 0 is 
equivalent to (Bl i, thus providing an equivalent form of the Peres-Horodecki-Simon criterion [24] of entanglement. ■ 
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C. Several Auxiliary Facts on Hermitian Matrices 

For ease of reference, we mention the following properties of Hermitian matrices S 
i) S 0 if and only if 

D a C C 

0 0 . 


Re.S’ + dm.S': 


ReS 

ImS 


—ImS 
ReS 


(Cl) 


ii) S tp 0 if and only if S 0 0. 

iii) Suppose S 0 0 is a (4 x 4)-matrix with ImS := h 0 ©, and 0 £ A 2 . Then all the principal minors of the matrix 

A := ReS + i 


'0 

0 ’ 


An 

A 12 

0 

-0 


A 12 

a 22 


(C2) 

are nonnegative, except probably det A. The latter determinant can be computed in terms of the (2 x 2)-blocks of A as 

det A = det A 11 det (A 2 2 - A^Af, A 1 2 ) 

= detA 22 det (An - A| 2 Aj 2 A} 2 ). (C3) 


Proof: (i) is verified by inspection. In view of (i), 
r t o’ 


ii) can be proved by left and right multiplying both sides of (Cl I 


. To prove (iii), note that the matrix A in (C2i shares with S or S all principal submatrices up to order 3. These 


by 

submatrices inherit positive semi-definiteness from S and S, and hence, the corresponding minors of A are nonnegative. 
Therefore, only the fourth order minor detA can be negative. Now, elementary matrix transformations lead to 


/ 

o' 

'All 

Ai 2 

'/ -Vi a i 2 ' 

- A 12 A lt 

I 

A 12 

A 22 

0 I 


Aj2 A} 2 AjjAii A 22 A| 2 A [' | A12 


All 

10 Ml 
An 0 

0 A 22 —AfjA+Aia 


A12 — An At, A12 


Here, use is made of the properties that A^AnA^ = A M and A = A* together with 

(AnA+ — 7)A 12 = (SnS+ -I)S n = 0, 


(C4) 


(C5) 


where the latter follows from An = Sn, Ai 2 = Si 2 and the assumption that S A 0 [31], Indeed, every S 0 0 has a unique 
square root \/S 0 0 whose columns can be partitioned according to the blocks of S. More precisely, let \/S := \L\ L 2 ], 

and hence. 


S = 


Sll 

S\2 


L*\L\ 

l\l 2 ' 

Sn 

S 22 


L%L X 

l* 2 l 2 


(C6) 


Consider the polar decomposition of L\ = <$ >5, where <f> has orthonormal columns and 3 = E* A 0. Then S\ 1 = L\L\ = E*S = 
E 2 and Si 2 = L\Li = E<t>*L 2 in view of (C 61 . This implies 5 h5j 1 5'i 2 = E 2 (E 2 ) + E<t>*L 2 = EE + S<I>*L 2 = S\ 2 , thus establishing 
(C5 1 . Here, use is also made of E + = E*(EE*) + , the fact that E is Hermitian and the definition of the generalized inverse. 


Now, by evaluating the determinant of both sides of (C4i, it follows that detA is indeed representable by the first of the 
equalities (|C3|>. The second of them is established in a similar fashion. ■ 
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